function [irf,irfse] = estimate_lp_iv(depvar,indepvar,shock)

irf = NaN(17,1);
irfse = NaN(17,1);

shock = shock./std(shock(~isnan(shock)));
x = indepvar-lagmatrix(indepvar,1);

for h=0:16
    y = lagmatrix(depvar,-h)-lagmatrix(depvar,1);
    X = [x     lagmatrix(depvar,1)-lagmatrix(depvar,2) ones(size(x))];
    Z = [shock lagmatrix(depvar,1)-lagmatrix(depvar,2) ones(size(shock))];
    [beta,Var] = nwest_iv(y,X,Z,h+1);
    irf(h+1)=beta(1);
    irfse(h+1)=sqrt(Var(1,1));
end

end
